#!/usr/bin/python
# GM_LumFuncBochanski.py
'''Create realistic populations of stars, up to a constant.'''

import numpy

VariableNames = ['z', 'R', 'r', 'R0', 'BoR0thin', 'BoR0thick', 'Bof', 'galactocentric', 'scale', 'LumFuncList', 'delMvList']

def lumfunc(simpars, phypars):
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    for i in BoCorrectedValues:
        cmd = "%s = BoCorrectedValues['%s']" % (i,i)
        exec cmd
    # R dependence of density
    RCorrThin = numpy.exp( -(R - R0)/BoR0thin )
    RCorrThick = numpy.exp( -(R - R0)/BoR0thick )
    # z dependence of density
    zCorrThin = numpy.exp( -z/BoZ0thin )
    zCorrThick = numpy.exp( -z/BoZ0thick )
    # Combine
    ThinDisk = numpy.multiply(RCorrThin, zCorrThin[:,numpy.newaxis,:])
    ThickDisk = numpy.multiply(RCorrThick, zCorrThick[:,numpy.newaxis,:])
    TotalDisk = Bof*ThinDisk + (1 - Bof)*ThickDisk
    # Optional cutoff at r>scale to prevent assymetry about galactic center
    if galactocentric == True:
        rmask = numpy.greater(scale, r)
        TotalDisk = TotalDisk*rmask
        del rmask
    elif galactocentric == False:
        pass
    else:
        print "You're a cucumber!  ...No seriously, check the value you entered for 'galactocentric'.  It should be boolean."
    distriblist = []
    del RCorrThin, RCorrThick, zCorrThin, zCorrThick, ThinDisk, ThickDisk
    for i in range(0, len(MvList)):
        rho0Mv = LumFuncList[i]
        binwidth = delMvList[i]
        thisdistrib = binwidth*rho0Mv*numpy.array(TotalDisk)
        distriblist.append(thisdistrib)
    del TotalDisk
    return distriblist
